
## read in primary data
dists <- readRDS("temp/shooting_demos.rds") |> 
  mutate(date = as.Date(date),
         id = id4, ## id (and distance, in next line) of nearest 4-victim shooting
         dist = dist4) |> 
  filter(!is.na(id),
         id %in% (readRDS("temp/geocoded_shootings_from_gva.rds") |> 
                    filter(home + dv == 0))$id ## drop if in the home or domestic violence
  )

bws <- readRDS("temp/primary_out_data.rds") |> 
  filter(p <= 10)

## loop over thresholds
run_rob <- function(threshold){
  library(data.table)
  library(tidyverse)
  library(ebal)
  library(rdrobust)
  ##keep the observations within the threshold closest to election day
  set_pre <- dists[dists$dist <= threshold & dists$pre,
                   .SD[date %in% max(date)], by=list(year, GEOID)]
  set_pre <- set_pre[, .SD[1], by = .(year, GEOID)] %>% 
    mutate(treated = T)
  
  ##keep the observations within the threshold closest to election day
  set_post <- dists[dists$dist <= threshold & !dists$pre,
                    .SD[date %in% min(date)], by=list(year, GEOID)]
  set_post <- set_post[!(paste0(set_post$GEOID, set_post$year) %in% paste0(set_pre$GEOID, set_pre$year)),
                       .SD[1], by = .(year, GEOID)] %>% 
    mutate(treated = F)
  
  full_treat <- bind_rows(set_pre, set_post) %>% 
    select(GEOID, id, date, dist, year, turnout,
           median_income, nh_white, nh_black, median_age, pop_dens, latino, asian,
           some_college, turnout_pre, treated) %>% 
    mutate(d2 = ifelse(year == "2020", as.integer(date - as.Date("2020-11-03")),
                       as.integer(date - as.Date("2016-11-08"))),
           t16 = year == "2016")
  
  full_treat <- full_treat %>%
    mutate(change_to = turnout - turnout_pre)
  ########################################
  
  ## main model with narrowest bandwidths
  la <- rdrobust(y = full_treat$turnout, x = full_treat$d2, p = 1, c = 0, cluster = full_treat$id,
                 covs = select(full_treat,
                               latino, nh_white, asian,
                               nh_black, median_income, median_age,
                               pop_dens, turnout_pre, t16, some_college), h = min(bws$bw))
  ## main model with widest bandwidths
  lb <- rdrobust(y = full_treat$turnout, x = full_treat$d2, p = 1, c = 0, cluster = full_treat$id,
                 covs = select(full_treat,
                               latino, nh_white, asian,
                               nh_black, median_income, median_age,
                               pop_dens, turnout_pre, t16, some_college), h = max(bws$bw))

  
  ## rdit on first difference in turnout (note the changed dependent variable, different weight, and
  ## noninclusion of turnout_pre in the model)
  l4 <- rdrobust(y = full_treat$change_to, x = full_treat$d2, p = 1, c = 0, cluster = full_treat$id,
                 covs = select(full_treat,
                               latino, nh_white, asian,
                               nh_black, median_income, median_age,
                               pop_dens, t16, some_college))
  
  l4b <- rdrobust(y = full_treat$turnout_pre, x = full_treat$d2, p = 1, c = 0, cluster = full_treat$id,
                  covs = select(full_treat,
                                latino, nh_white, asian,
                                nh_black, median_income, median_age,
                                pop_dens, t16, some_college))
  
  ## grab nonparametric bandwidth for alternate specifications
  bw <- rdbwselect(y = full_treat$turnout, x = full_treat$d2, p = 1, c = 0, cluster = full_treat$id,
                   # weights = full_treat$weight,
                   covs = select(full_treat,
                                 latino, nh_white, asian,
                                 nh_black, median_income, median_age,
                                 pop_dens, turnout_pre, t16, some_college))
  
  ## rdit using double the nonparametric bandwidth
  l5 <- rdrobust(y = full_treat$turnout, x = full_treat$d2, p = 1, c = 0, cluster = full_treat$id,
                 covs = select(full_treat,
                               latino, nh_white, asian,
                               nh_black, median_income, median_age,
                               pop_dens, turnout_pre, t16, some_college), h = bw[["bws"]][1]*2)
  ## rdit using half the nonparametric bandwidth
  l5b <- rdrobust(y = full_treat$turnout, x = full_treat$d2, p = 1, c = 0, cluster = full_treat$id,
                 covs = select(full_treat,
                               latino, nh_white, asian,
                               nh_black, median_income, median_age,
                               pop_dens, turnout_pre, t16, some_college), h = bw[["bws"]][1]/2)
  ## rdit using nonparametric bandwidth on treated side, 60 days on untreated side
  l6 <- rdrobust(y = full_treat$turnout, x = full_treat$d2, p = 1, c = 0, cluster = full_treat$id,
                 covs = select(full_treat,
                               latino, nh_white, asian,
                               nh_black, median_income, median_age,
                               pop_dens, turnout_pre, t16, some_college), h = c(bw[["bws"]][1], 60))
  ## rdit using nonparametric bandwidth on treated side, 90 days on untreated side
  l7 <- rdrobust(y = full_treat$turnout, x = full_treat$d2, p = 1, c = 0, cluster = full_treat$id,
                 covs = select(full_treat,
                               latino, nh_white, asian,
                               nh_black, median_income, median_age,
                               pop_dens, turnout_pre, t16, some_college), h = c(bw[["bws"]][1], 90))
  ## rdit using nonparametric bandwidth on treated side, 180 days on untreated side
  l7b <- rdrobust(y = full_treat$turnout, x = full_treat$d2, p = 1, c = 0, cluster = full_treat$id,
                 covs = select(full_treat,
                               latino, nh_white, asian,
                               nh_black, median_income, median_age,
                               pop_dens, turnout_pre, t16, some_college), h = c(bw[["bws"]][1], 180))
  ## keep only the 2016 data
  o16 <- filter(full_treat, year == "2016")
  ## rdit on 2016 data alone
  l8 <- rdrobust(y = o16$turnout, x = o16$d2, p = 1, c = 0, cluster = o16$id,
                 covs = select(o16,
                               latino, nh_white, asian,
                               nh_black, median_income, median_age,
                               pop_dens, turnout_pre, some_college))
  ##keep only the 2020 data
  o20 <- filter(full_treat, year != "2016")
  #rdit on 2020 alone
  l9 <- rdrobust(y = o20$turnout, x = o20$d2, p = 1, c = 0, cluster = o20$id,
                 covs = select(o20,
                               latino, nh_white, asian,
                               nh_black, median_income, median_age,
                               pop_dens, turnout_pre, some_college))
  
  no_gangs <- filter(full_treat,
                     id %in% (readRDS("temp/geocoded_shootings_from_gva.rds") |> 
                                filter(gang == 0))$id ## drop if gang
  )
  #rdit on 2020 alone
  l10 <- rdrobust(y = no_gangs$turnout, x = no_gangs$d2, p = 1, c = 0, cluster = no_gangs$id,
                 covs = select(no_gangs,
                               latino, nh_white, asian,
                               nh_black, median_income, median_age,
                               pop_dens, turnout_pre, some_college))
  
  ## combine results of all these models
  f <- rbind(
    tibble(coef = la$coef,
           n = la$N_h[1] + la$N_h[2],
           bwidth = floor(la[["bws"]][2]),
           se = la$se, 
           pv = la$pv,
           p = threshold,
           t = "Narrowest",
           u = la[["ci"]][,2],
           l = la[["ci"]][,1]),
    tibble(coef = lb$coef,
           n = lb$N_h[1] + lb$N_h[2],
           bwidth = floor(lb[["bws"]][2]),
           se = lb$se, 
           pv = lb$pv,
           p = threshold,
           t = "Widest",
           u = lb[["ci"]][,2],
           l = lb[["ci"]][,1]),
    tibble(coef = l4$coef,
           n = l4$N_h[1] + l4$N_h[2],
           bwidth = floor(l4[["bws"]][2]),
           se = l4$se, 
           pv = l4$pv,
           p = threshold,
           t = "First Difference in Turnout",
           u = l4[["ci"]][,2],
           l = l4[["ci"]][,1]),
    tibble(coef = l4b$coef,
           n = l4b$N_h[1] + l4b$N_h[2],
           bwidth = floor(l4b[["bws"]][2]),
           se = l4b$se, 
           pv = l4b$pv,
           p = threshold,
           t = "Prior Turnout",
           u = l4b[["ci"]][,2],
           l = l4b[["ci"]][,1]),
    tibble(coef = l5$coef,
           n = l5$N_h[1] + l5$N_h[2],
           bwidth = floor(l5[["bws"]][2]),
           se = l5$se, 
           pv = l5$pv,
           p = threshold,
           t = "Double Bandwidth",
           u = l5[["ci"]][,2],
           l = l5[["ci"]][,1]),
    tibble(coef = l5b$coef,
           n = l5b$N_h[1] + l5b$N_h[2],
           bwidth = floor(l5b[["bws"]][2]),
           se = l5b$se, 
           pv = l5b$pv,
           p = threshold,
           t = "Half Bandwidth",
           u = l5b[["ci"]][,2],
           l = l5b[["ci"]][,1]),
    tibble(coef = l6$coef,
           n = l6$N_h[1] + l6$N_h[2],
           se = l6$se, 
           pv = l6$pv,
           p = threshold,
           t = "Nonpara, 60",
           bwidth = paste0("(", floor(l6[["bws"]][1]), ", ", floor(l6[["bws"]][3]), ")"),
           u = l6[["ci"]][,2],
           l = l6[["ci"]][,1]),
    tibble(coef = l7$coef,
           n = l7$N_h[1] + l7$N_h[2],
           bwidth = paste0("(", floor(l7[["bws"]][1]), ", ", floor(l7[["bws"]][3]), ")"),
           se = l7$se, 
           pv = l7$pv,
           p = threshold,
           t = "Nonpara, 90",
           u = l7[["ci"]][,2],
           l = l7[["ci"]][,1]),
    tibble(coef = l7b$coef,
           n = l7b$N_h[1] + l7b$N_h[2],
           bwidth = paste0("(", floor(l7b[["bws"]][1]), ", ", floor(l7b[["bws"]][3]), ")"),
           se = l7b$se, 
           pv = l7b$pv,
           p = threshold,
           t = "Nonpara, 180",
           u = l7b[["ci"]][,2],
           l = l7b[["ci"]][,1]),
    tibble(coef = l8$coef,
           n = l8$N_h[1] + l8$N_h[2],
           bwidth = floor(l8[["bws"]][2]),
           se = l8$se, 
           pv = l8$pv,
           p = threshold,
           t = "Only 2016",
           u = l8[["ci"]][,2],
           l = l8[["ci"]][,1]),
    tibble(coef = l9$coef,
           n = l9$N_h[1] + l9$N_h[2],
           bwidth = floor(l9[["bws"]][2]),
           se = l9$se, 
           pv = l9$pv,
           p = threshold,
           t = "Only 2020",
           u = l9[["ci"]][,2],
           l = l9[["ci"]][,1]),
    tibble(coef = l10$coef,
           n = l10$N_h[1] + l10$N_h[2],
           bwidth = floor(l10[["bws"]][2]),
           se = l10$se, 
           pv = l10$pv,
           p = threshold,
           t = "No Gangs",
           u = l10[["ci"]][,2],
           l = l10[["ci"]][,1])
  )
}

cl <- makeCluster(12)  
registerDoParallel(cl)

clusterExport(cl, list("run_rob", "dists", "bws"))

out <- rbindlist(parLapply(cl, unique(bws$p), fun = run_rob))

saveRDS(out, "temp/alt_rdds2.rds")

########### figure s14
out <- readRDS("temp/alt_rdds2.rds") %>% 
  filter(t == "First Difference in Turnout")

out$estimate <- rep(c('Traditional','Bias-Adjusted','Robust'),nrow(out)/3)

out <- mutate_at(out, vars(coef, l, u), ~. * -1)

different_dists <- ggplot(filter(out, estimate == "Robust"),
                          aes(x = p, y = coef, ymin = l, ymax = u)) +
  geom_point() +
  geom_errorbar(width = 0) + 
  theme_bc(base_family = "LM Roman 10", base_size = 10,
           plot.margin = margin()) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_x_continuous(breaks = seq(0, 10, 2)) +
  labs(y = "Local average treatment effect", x = "Radius around shooting (miles)")
different_dists

fwrite(out |> 
         filter(estimate == "Robust") |> 
         select(distance = p,
                estimate = coef,
                lower_bound = l,
                upper_bound = u), "data_for_figures/figure_A8.csv")

saveRDS(different_dists, "temp/first_diff_plot.rds")
for(i in outpaths){
  ggsave(paste0(i, "first_difference_plot.png"), height = 4, width = 3, units = "in")
}
######################
#create figure s11
out <- readRDS("temp/alt_rdds2.rds") %>% 
  filter(t %in% c('Half Bandwidth', 'Double Bandwidth','Nonpara, 60','Nonpara, 90', 'Nonpara, 180')) %>% 
  mutate(t = ifelse(t == "Nonpara, 60", "60 Days",
                    ifelse(t == "Nonpara, 90", "90 Days",
                           ifelse(t == "Nonpara, 180", "180 Days", gsub(" ", "\n", t)))))

out$estimate <- rep(c('Traditional','Bias-Adjusted','Robust'),nrow(out)/3)

out <- mutate_at(out, vars(coef, l, u), ~. * -1)

out$estimate <- factor(out$estimate, levels = c('Traditional','Bias-Adjusted','Robust'))
out$t <- factor(out$t, levels = c('Half\nBandwidth', 'Double\nBandwidth','60 Days','90 Days', '180 Days'))

different_dists <- ggplot(filter(out, estimate == "Robust"),
                          aes(x = p, y = coef, ymin = l, ymax = u)) +
  facet_grid(t~.) +
  geom_point() +
  geom_errorbar(width = 0) + 
  theme_bc(base_family = "LM Roman 10", base_size = 10,
           plot.margin = margin()) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_x_continuous(breaks = seq(0, 10, 2)) +
  labs(y = "Local average treatment effect", x = "Radius around shooting (miles)")
different_dists

fwrite(out |> 
         filter(estimate == "Robust") |> 
         select(bandwidth = t,
                distance = p,
                estimate = coef,
                lower_bound = l,
                upper_bound = u), "data_for_figures/figure_A5.csv")

saveRDS(different_dists, "temp/alt_bws_rdd.rds")
for(i in outpaths){
  ggsave(paste0(i, "different_bandwidths.png"), height = 7.5, width = 3, units = "in")
}

######################
## create figure s6
out <- readRDS("temp/alt_rdds2.rds") %>% 
  filter(t %in% c('Only 2016', 'Only 2020'))

out$estimate <- rep(c('Traditional','Bias-Adjusted','Robust'),nrow(out)/3)

out <- mutate_at(out, vars(coef, l, u), ~. * -1)

out$estimate <- factor(out$estimate, levels = c('Traditional','Bias-Adjusted','Robust'))
# out$t <- factor(out$t, levels = c('Double Bandwidth','30 Days','60 Days'))

different_dists <- ggplot(out,
                          aes(x = p, y = coef, ymin = l, ymax = u)) +
  facet_grid(t~estimate) +
  geom_point() +
  geom_errorbar(width = 0) + 
  theme_bc(base_family = "LM Roman 10", base_size = 13) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_x_continuous(breaks = seq(0, 10, 2)) +
  labs(y = "Local average treatment effect", x = "Radius around shooting (miles)")
different_dists

saveRDS(different_dists, "temp/individual_years.rds")
for(i in outpaths){
  ggsave(paste0(i, "individual_years.png"), height = 4, width = 6, units = "in")
}

######################
## create figure s13
out <- readRDS("temp/alt_rdds2.rds") %>% 
  filter(t %in% c("Prior Turnout"))

out$estimate <- rep(c('Traditional','Bias-Adjusted','Robust'),nrow(out)/3)

out <- mutate_at(out, vars(coef, l, u), ~. * -1)

out$estimate <- factor(out$estimate, levels = c('Traditional','Bias-Adjusted','Robust'))
# out$t <- factor(out$t, levels = c('Double Bandwidth','30 Days','60 Days'))

different_dists <- ggplot(filter(out, estimate == "Robust"),
                          aes(x = p, y = coef, ymin = l, ymax = u)) +
  geom_point() +
  geom_errorbar(width = 0) + 
  theme_bc(base_family = "LM Roman 10", base_size = 10,
           plot.margin = margin()) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_x_continuous(breaks = seq(0, 10, 2)) +
  labs(y = "Local average treatment effect", x = "Radius around shooting (miles)")
different_dists
saveRDS(different_dists, "temp/placebo_prior.rds")

fwrite(out |> 
         filter(estimate == "Robust") |> 
         select(distance = p,
                estimate = coef,
                lower_bound = l,
                upper_bound = u), "data_for_figures/figure_A7.csv")

for(i in outpaths){
  ggsave(paste0(i, "placebo_prior_turnout.png"), height = 4, width = 3, units = "in")
}
######################
## create figure s12
out <- readRDS("temp/alt_rdds2.rds") %>% 
  filter(t %in% c("Narrowest", "Widest"))

out$estimate <- rep(c('Traditional','Bias-Adjusted','Robust'),nrow(out)/3)

out <- mutate_at(out, vars(coef, l, u), ~. * -1)

out$estimate <- factor(out$estimate, levels = c('Traditional','Bias-Adjusted','Robust'))
# out$t <- factor(out$t, levels = c('Double Bandwidth','30 Days','60 Days'))

different_dists <- ggplot(filter(out, estimate == "Robust"),
                          aes(x = p, y = coef, ymin = l, ymax = u)) +
  facet_grid(t~.) +
  geom_point() +
  geom_errorbar(width = 0) + 
  theme_bc(base_family = "LM Roman 10", base_size = 10,
           plot.margin = margin()) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_x_continuous(breaks = seq(0, 10, 2)) +
  labs(y = "Local average treatment effect", x = "Radius around shooting (miles)")
different_dists

fwrite(out |> 
         filter(estimate == "Robust") |> 
         select(bandwidth = t,
                distance = p,
                estimate = coef,
                lower_bound = l,
                upper_bound = u), "data_for_figures/figure_A6.csv")

saveRDS(different_dists, "temp/same_bws.rds")
for(i in outpaths){
  ggsave(paste0(i, "same_bandwidth_throughout.png"), height = 4, width = 3, units = "in")
}

###############################
###############################
###############################
###############################
###############################


dists <- readRDS("temp/shooting_demos.rds") |> 
  mutate(date = as.Date(date),
         id = id4, ## id (and distance, in next line) of nearest 4-victim shooting
         dist = dist4) |> 
  filter(!is.na(id),
         id %in% (readRDS("temp/geocoded_shootings_from_gva.rds") |> 
                    filter(home + dv + gang == 0))$id ## drop if in the home or domestic violence or gang
  )

### loop over thresholds for RDITS
run_rdd <- function(threshold){
  library(data.table)
  library(tidyverse)
  library(rdrobust)
  print(threshold)
  ##keep the observations within the threshold closest to election day
  ## using data.table notation to speed this up
  set_pre <- dists[dists$dist <= threshold & dists$pre,
                   .SD[date %in% max(date)], by=list(year, GEOID)]
  set_pre <- set_pre[, .SD[1], by = .(year, GEOID)] %>% 
    mutate(treated = T)
  
  ##keep the observations within the threshold closest to election day
  set_post <- dists[dists$dist <= threshold & !dists$pre,
                    .SD[date %in% min(date)], by=list(year, GEOID)]
  set_post <- set_post[!(paste0(set_post$GEOID, set_post$year) %in% paste0(set_pre$GEOID, set_pre$year)),
                       .SD[1], by = .(year, GEOID)] %>% 
    mutate(treated = F)
  
  full_treat <- bind_rows(set_pre, set_post) %>% 
    select(GEOID, id, date, dist, year, turnout,
           median_income, nh_white, nh_black, median_age, pop_dens, latino, asian,
           some_college, turnout_pre, treated) %>% 
    mutate(d2 = ifelse(year == "2020", as.integer(date - as.Date("2020-11-03")),
                       as.integer(date - as.Date("2016-11-08"))),
           t16 = year == "2016")
  

  ########################################
  
  l <- rdrobust(y = full_treat$turnout, x = full_treat$d2, p = 1, c = 0, cluster = full_treat$id,
                covs = select(full_treat,
                              latino, nh_white, asian,
                              nh_black, median_income, median_age,
                              pop_dens, turnout_pre, t16, some_college))
  
  f <- tibble(coef = l$coef,
              se = l$se, 
              pv = l$pv,
              p = threshold,
              n = l$N_h[1] + l$N_h[2],
              bw = l[["bws"]][1],
              u = l[["ci"]][,2],
              l = l[["ci"]][,1])
}

cl <- makeCluster(12)  
registerDoParallel(cl)

clusterExport(cl, list("run_rdd", "dists"))
out <- rbindlist(parLapply(cl, c(0.25, 0.5, 0.75, seq(1, 10, 1)), fun = run_rdd))

saveRDS(out, "temp/no_gangs_out_data.rds")
## create figure s15
out <- readRDS("temp/no_gangs_out_data.rds")
out$estimate <- rep(c('Traditional','Bias-Adjusted','Robust'),nrow(out)/3)

out <- mutate_at(out, vars(coef, l, u), ~. * -1)

out$estimate <- factor(out$estimate, levels = c('Traditional','Bias-Adjusted','Robust'))

different_dists <- ggplot(filter(out, estimate == "Robust"),
                          aes(x = p, y = coef, ymin = l, ymax = u)) +
  geom_point() +
  geom_errorbar(width = 0) + 
  theme_bc(base_family = "LM Roman 10", base_size = 10,
           plot.margin = margin()) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_x_continuous(breaks = seq(0, 10, 2)) +
  labs(y = "Local average treatment effect", x = "Radius around shooting (miles)")
different_dists

fwrite(out |> 
         filter(estimate == "Robust") |> 
         select(distance = p,
                estimate = coef,
                lower_bound = l,
                upper_bound = u), "data_for_figures/figure_A9.csv")

saveRDS(different_dists, "temp/different_dists_no_gangs.rds")

for(i in outpaths){
  ggsave(paste0(i, "no_gangs_results.png"), height = 4, width = 3, units = "in")
}

###############################
###############################
###############################
###############################
###############################


dists <- readRDS("temp/shooting_demos.rds") |> 
  mutate(date = as.Date(date),
         id = id4, ## id (and distance, in next line) of nearest 4-victim shooting
         dist = dist4) |> 
  filter(!is.na(id),
         id %in% (readRDS("temp/geocoded_shootings_from_gva.rds") |> 
                    filter(home + dv + hate + political + sex == 0))$id ## drop if in the home or domestic violence or gang
  )



### loop over thresholds for RDITS
run_rdd <- function(threshold){
  library(data.table)
  library(tidyverse)
  library(rdrobust)
  print(threshold)
  ##keep the observations within the threshold closest to election day
  ## using data.table notation to speed this up
  set_pre <- dists[dists$dist <= threshold & dists$pre,
                   .SD[date %in% max(date)], by=list(year, GEOID)]
  set_pre <- set_pre[, .SD[1], by = .(year, GEOID)] %>% 
    mutate(treated = T)
  
  ##keep the observations within the threshold closest to election day
  set_post <- dists[dists$dist <= threshold & !dists$pre,
                    .SD[date %in% min(date)], by=list(year, GEOID)]
  set_post <- set_post[!(paste0(set_post$GEOID, set_post$year) %in% paste0(set_pre$GEOID, set_pre$year)),
                       .SD[1], by = .(year, GEOID)] %>% 
    mutate(treated = F)
  
  full_treat <- bind_rows(set_pre, set_post) %>% 
    select(GEOID, id, date, dist, year, turnout,
           median_income, nh_white, nh_black, median_age, pop_dens, latino, asian,
           some_college, turnout_pre, treated) %>% 
    mutate(d2 = ifelse(year == "2020", as.integer(date - as.Date("2020-11-03")),
                       as.integer(date - as.Date("2016-11-08"))),
           t16 = year == "2016")
  
  
  ########################################
  
  l <- rdrobust(y = full_treat$turnout, x = full_treat$d2, p = 1, c = 0, cluster = full_treat$id,
                covs = select(full_treat,
                              latino, nh_white, asian,
                              nh_black, median_income, median_age,
                              pop_dens, turnout_pre, t16, some_college))
  
  f <- tibble(coef = l$coef,
              se = l$se, 
              pv = l$pv,
              p = threshold,
              n = l$N_h[1] + l$N_h[2],
              bw = l[["bws"]][1],
              u = l[["ci"]][,2],
              l = l[["ci"]][,1])
}

cl <- makeCluster(12)  
registerDoParallel(cl)

clusterExport(cl, list("run_rdd", "dists"))
out <- rbindlist(parLapply(cl, c(0.25, 0.5, 0.75, seq(1, 10, 1)), fun = run_rdd))

saveRDS(out, "temp/no_political_hate_sex_out_data.rds")

out <- readRDS("temp/no_political_hate_sex_out_data.rds")
out$estimate <- rep(c('Traditional','Bias-Adjusted','Robust'),nrow(out)/3)

out <- mutate_at(out, vars(coef, l, u), ~. * -1)

out$estimate <- factor(out$estimate, levels = c('Traditional','Bias-Adjusted','Robust'))
## figure s19
different_dists <- ggplot(filter(out, estimate == "Robust"),
                          aes(x = p, y = coef, ymin = l, ymax = u)) +
  # facet_grid(~estimate) +
  geom_point() +
  geom_errorbar(width = 0) + 
  theme_bc(base_family = "LM Roman 10", base_size = 10,
           plot.margin = margin()) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_x_continuous(breaks = seq(0, 10, 2)) +
  labs(y = "Local average treatment effect", x = "Radius around shooting (miles)")
different_dists

fwrite(out |> 
         filter(estimate == "Robust") |> 
         select(distance = p,
                estimate = coef,
                lower_bound = l,
                upper_bound = u), "data_for_figures/figure_no_political.csv")

saveRDS(different_dists, "temp/different_dists_no_political.rds")

for(i in outpaths){
  ggsave(paste0(i, "no_political_results.png"), height = 4, width = 3, units = "in")
}
